Marxan in R

Tutorial for ESM 270 Conservation Planning

Introduction

Designing reserve networks requires considering both ecological and socioeconomic factors. Marxan is a commonly used spatial planning tool that quantifies these factors to generate efficient reserve networks. First we will run through a simple example analysis, and then we will use species and parcel data from a 2014-15 Bren Group Project in the Morro Bay Watershed to examine the optimate regional reserve network.

Priotitizr

Prioritizr is an R package that, like Marxan, can guide systematic reserve design by solving conservation prioritization problems. In this tutorial, we will be using the marxan_problem() function which can read Marxan input data. More information on the prioritizr package can be found at https://mran.microsoft.com/snapshot/2018-03-04/web/packages/prioritizr/vignettes/prioritizr_basics.html.

A really useful tutorial in using prioritizr, specifically in using summed solutions, can be found here: https://prioritizr.net/articles/tasmania.html. This resource helped frame this tutorial.

Gurobi

For the solve() function, you must have a solver installed. For the the process with summing and multiple runs process, this must be gurobi. This requires following the instructions below.

1. Download gurobi (https://www.gurobi.com/downloads/gurobi-optimizer-eula/)
2. Request an academic license (https://www.gurobi.com/downloads/end-user-license-agreement-academic/)
3. Verify license in terminal by copying and pasting information in your account > your licenses
4. Must install r package withinstall.packages('/Library/gurobi902/mac64/R/gurobi_9.0-2_R_3.6.1.tgz', repos=NULL)
5. must install slam package with `install.packages(“slam”)
6. add gurobi and slam libraries

Getting Started

Attach Packages

library(here)
library(janitor)
library(prioritizr) # marxan
library(sf) # spatial features
library(slam) # to use gurobi
library(gurobi) # solver
library(ggmap) # basemaps
library(patchwork) 
library(tidyverse)

If you do not have any of these packages installed, run the following in the console: install.packages("package name")

Note: this will not work for installing gurobi. Instructions are in the introduction section.

Example Analysis

Familiarize yourself with the format of Marxan input files (all are text files that can be opened with a text editor) using this small example.

Data

• spec.dat: species file that contains 24 elements with varying targets expressed as the number of sites protected, constant species penalty factor, and a minimum of 2 occurrences.

• pudata.dat: planning unit file that lists 99 sites with variable site costs.

• PUvCVr.dat: planning unit vs. species file that lists the amount of each species in each planning unit in a codensed format.

• bound.dat: boundary length file that provides the ‘effective length’ of shared boundary between pairs of planning units.

Read in Data

Read in data and select necessary columns for marxan_problem(). Find more information with ?marxan_problem().

spec_samp <- read_tsv(here("sample-data", "spec.dat")) %>% 
  rename("amount" = "target")
# with prioritizr, target is called "prop" (relative) or "amount" (absolute)

pu_samp <- read_tsv(here("sample-data","pudata.dat"))

puvsp_samp <- read_tsv(here("sample-data", "PUvCVr.dat"))

bound_samp <- read_tsv(here("sample-data", "bound.dat"))

Run Marxan

The function marxan_problem() in the prioritizr package gives us the “canned” approach that works for our purposes. If more customizations are desired, feel free to explore the problem() function instead.

Arguments for marxan_problem (for more information, see ?marxan_problem):

• x: planning units to use in reserve design and their cost (x = pu_samp)
• spec: conservation features. Must contain columns for id, name, and prop/amount (spec = spec_samp)
• puvspr: the amount of each feature in each planning unit. Must contain columns for pu, species, and amount. (puvspr = puvsp_samp)
• bound: boundary data (bound = bound_samp)
• blm: boundary length modifier. (blm = 0.2)

Also note that in solver, we will use 100 runs (number_solutions = 100)

marx_samp <- marxan_problem(x = pu_samp,
                              spec = spec_samp,
                              puvspr = puvsp_samp,
                              bound = bound_samp,
                              blm = 0.2)

marx_prob_samp <- marx_samp %>% 
  add_gurobi_solver(gap = 0.15) %>% 
  add_pool_portfolio(method = 2,
                     number_solutions = 10) %>% 
  add_absolute_targets(2)

# gap = optimality gap. 0.15 gives us results closest to that of using marxan software. Learn more under ?add_gurobi_solver 
# add_absolute_targets sets targets for the actual value of features in the study area that need to be represented in the prioritization. Learn more under ?add_absolute_targets. Other option is add_relative_targets(). 

marx_samp_soln <- solve(marx_prob_samp)

## Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (mac64)
## Optimize a model with 380 rows, 277 columns and 2116 nonzeros
## Model fingerprint: 0x3b987834
## Variable types: 0 continuous, 277 integer (277 binary)
## Coefficient statistics:
##   Matrix range     [1e+00, 8e+01]
##   Objective range  [4e-01, 2e+01]
##   Bounds range     [1e+00, 1e+00]
##   RHS range        [2e+00, 2e+00]
## Found heuristic solution: objective 1457.2100000
## Found heuristic solution: objective 130.1700000
## Presolve removed 1 rows and 0 columns
## Presolve time: 0.01s
## Presolved: 379 rows, 277 columns, 2038 nonzeros
## Variable types: 0 continuous, 277 integer (277 binary)
## Presolve removed 6 rows and 0 columns
## Presolved: 373 rows, 277 columns, 1475 nonzeros
## 
## 
## Root relaxation: objective 7.421647e+01, 29 iterations, 0.00 seconds
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
##      0     0   74.21647    0    7  130.17000   74.21647  43.0%     -    0s
## H    0     0                     124.5600000   74.21647  40.4%     -    0s
## H    0     0                     123.5000000   74.21647  39.9%     -    0s
## H    0     0                      89.3300000   74.21647  16.9%     -    0s
## H    0     0                      83.3400000   74.21647  10.9%     -    0s
## H    0     0                      79.4600000   74.21647  6.60%     -    0s
##      0     0   76.42240    0   11   79.46000   76.42240  3.82%     -    0s
##      0     0   77.15375    0   11   79.46000   77.15375  2.90%     -    0s
##      0     0   77.23600    0   10   79.46000   77.23600  2.80%     -    0s
##      0     0   77.26632    0   12   79.46000   77.26632  2.76%     -    0s
##      0     0   77.41333    0   10   79.46000   77.41333  2.58%     -    0s
##      0     0   78.41695    0   12   79.46000   78.41695  1.31%     -    0s
##      0     0   78.55519    0   12   79.46000   78.55519  1.14%     -    0s
##      0     0   78.72242    0   11   79.46000   78.72242  0.93%     -    0s
##      0     0   79.05773    0   14   79.46000   79.05773  0.51%     -    0s
##      0     0   79.09317    0   13   79.46000   79.09317  0.46%     -    0s
##      0     0   79.10461    0   13   79.46000   79.10461  0.45%     -    0s
##      0     0   79.11786    0   14   79.46000   79.11786  0.43%     -    0s
##      0     0   79.11878    0   14   79.46000   79.11878  0.43%     -    0s
##      0     0   79.25167    0   13   79.46000   79.25167  0.26%     -    0s
##      0     0   79.34043    0   13   79.46000   79.34043  0.15%     -    0s
##      0     0          -    0        79.46000   79.46000  0.00%     -    0s
## 
## Optimal solution found at node 0 - now completing solution pool...
## 
##     Nodes    |    Current Node    |      Pool Obj. Bounds     |     Work
##              |                    |   Worst                   |
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
##      0     0          -    0       123.50000   79.46000  35.7%     -    0s
##      0     0          -    0       123.50000   79.46000  35.7%     -    0s
##      0     2          -    0       123.50000   79.46000  35.7%     -    0s
## 
## Cutting planes:
##   Gomory: 2
##   Cover: 2
##   MIR: 3
##   StrongCG: 1
##   Zero half: 1
## 
## Explored 26 nodes (258 simplex iterations) in 0.07 seconds
## Thread count was 1 (of 4 available processors)
## 
## Solution count 10: 79.46 81.29 83.34 ... 89.61
## 
## Optimal solution found (tolerance 1.50e-01)
## Best objective 7.946000000000e+01, best bound 7.946000000000e+01, gap 0.0000%

Morro Bay Watershed

We will be using a case study (and the data they prepared) from a past group project on the Morro Bay National Estuary Program.

Data

Species data

• MorroBay_species_polyons: shapefiles of species locations
• MorroBay_species_pts: point data of species locations
• MorroBay_habitats: shapefile of location of key habitats
• morro-bay-spec.csv: information on species/habitat ID (the conservation feature), target value for that species (these values are 30% of the # of parcels containing that species or 10% for habitats), and the species name. The ‘spf’ column is the ‘species penalty factor’ and is a multiplier that determines the penalty if the conservation objective is not met
• spec-name-status.csv: information on species status

Parcel data

• morro-bay-pu.csv: includes information on the planning unit ID, the cost of that planning unit, and each parcel’s status (value 0 = available for prioritization, value 3 = locking in, 3 = locked out)
• morro-bay-puvspr.csv: includes information on the amount (presence/absence) of each species in each planning unit (pu)

Analyze All Species

Read in Data

Species and pu information

spec <- read_csv(here("morro-bay-data", "morro-bay-spec.csv")) %>% 
  head(140) %>%  # mine reads in extra blank rows - skip if yours does not 
  rename("amount" = "target") # target is called "prop" (relative) or "amount" (absolute)

pu <- read_csv(here("morro-bay-data", "morro-bay-pu.csv")) %>% 
  select(1:3) # mine reads in an extra blank column - skip if yours does not

puvsp <- read_csv(here("morro-bay-data", "morro-bay-puvspr.csv")) %>% 
  select(1:3) %>% 
  head(11849)

status <- read_csv(here("morro-bay-data", "spec-name-status.csv")) %>% 
  select(1:3) %>% 
  head(140)

Polygons

parcels <- read_sf(dsn = here("morro-bay-data"), layer = "MorroBay_parcels") %>% 
  clean_names()

ggplot(data = parcels) +
  geom_sf()

Run Marxan

The function marxan_problem() in the prioritizr package gives us the “canned” approach that works for our purposes. If more customizations are desired, feel free to explore the problem() function instead.

Arguments for marxan_problem (for more information, see ?marxan_problem):

• x: planning units to use in reserve design and their cost (x = pu)
• spec: conservation features. Must contain columns for id, name, and prop/amount (spec = spec)
• puvspr: the amount of each feature in each planning unit. Must contain columns for pu, species, and amount. (puvspr = puvsp)
• bound: boundary data (we will not use this; bound = NULL)
• blm: boundary length modifier. (blm = 0)

Also note that in solver, we will use 100 runs (number_solutions = 100)

marx_mb <- marxan_problem(x = pu,
                         spec = spec, 
                         puvspr = puvsp, 
                         bound = NULL,
                         blm = 0)

marx_mb_problem <- marx_mb %>% 
  add_gurobi_solver(gap = 0.15) %>% 
  add_pool_portfolio(method = 2, number_solutions = 100) 
# gap = optimality gap. 0.15 gives us results closest to that of using marxan software. Learn more under ?add_gurobi_solver 
# method 2 finds a specified number of solutions that are nearest to optimality. Learn more under ?add_pool_portfolio

marx_mb_soln <- solve(marx_mb_problem)

## Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (mac64)
## Optimize a model with 140 rows, 670 columns and 11849 nonzeros
## Model fingerprint: 0xd2c5bccb
## Variable types: 0 continuous, 670 integer (670 binary)
## Coefficient statistics:
##   Matrix range     [1e+00, 1e+00]
##   Objective range  [2e+00, 5e+06]
##   Bounds range     [1e+00, 1e+00]
##   RHS range        [3e-01, 2e+02]
## Found heuristic solution: objective 4.238988e+07
## Presolve removed 116 rows and 205 columns
## Presolve time: 0.02s
## Presolved: 24 rows, 465 columns, 3066 nonzeros
## Variable types: 0 continuous, 465 integer (465 binary)
## Presolve removed 1 rows and 0 columns
## Presolved: 23 rows, 465 columns, 3058 nonzeros
## 
## 
## Root relaxation: objective 2.840941e+07, 14 iterations, 0.00 seconds
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
## *    0     0               0    2.840941e+07 2.8409e+07  0.00%     -    0s
## 
## Optimal solution found at node 0 - now completing solution pool...
## 
##     Nodes    |    Current Node    |      Pool Obj. Bounds     |     Work
##              |                    |   Worst                   |
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
##      0     0          -    0               - 2.8409e+07      -     -    0s
##      0     0          -    0               - 2.8409e+07      -     -    0s
##      0     2          -    0               - 2.8409e+07      -     -    0s
## 
## Explored 1958 nodes (383 simplex iterations) in 0.18 seconds
## Thread count was 1 (of 4 available processors)
## 
## Solution count 100: 2.84094e+07 2.84097e+07 2.842e+07 ... 3.34128e+07
## 
## Optimal solution found (tolerance 1.50e-01)
## Best objective 2.840941400000e+07, best bound 2.840941400000e+07, gap 0.0000%

Sum Solutions

marx_mb_ssoln <- marx_mb_soln %>% 
  mutate(sum = rowSums(.[6:105])) %>% 
  select(id, cost, status, sum)

hist(marx_mb_ssoln$sum,
     main = "Selection Frequencies",
     xlab = "Number of runs that the unit was selected")

Visualize Results

Join polygons with output

ssoln_parcels <- inner_join(parcels, marx_mb_ssoln, by = "id")

ggplot(data = ssoln_parcels) +
  geom_sf(aes(fill = sum),
          color = "white",
          size = 0.05) + 
  scale_fill_gradient(low = "slategray2",
                      high = "navy") + 
  labs(fill = "Summed \nSolution") +
  theme_minimal()

Add basemap

This step is optional, but will make your map look better and is good to learn for further spatial analyses in R. Unfortunately, the options for basemaps in R are limited. Here, I have used the ggmaps package because it is the most accessible, however, it is still not very accessible compared to other R packages because you’ll need to get an API key. Instructions are here: https://cran.r-project.org/web/packages/ggmap/readme/README.html

For a ggmap cheatsheet, including the different basemaps in ggmaps, see: https://www.nceas.ucsb.edu/sites/default/files/2020-04/ggmapCheatsheet.pdf

morrobay <- get_map(location = c(lon = -120.7665, lat = 35.335),
                    zoom = 12,
                    maptype = "terrain-background", # background means no references, omit if want references
                    source = "google")

all_spec <- 
  ggmap(morrobay) +
    geom_sf(data = ssoln_parcels,
            aes(fill = sum),
          color = "white",
          size = 0.1,
          alpha = 0.85,
          inherit.aes = FALSE) + 
  coord_sf(crs = st_crs(4326)) +
  scale_fill_gradient(low = "slategray2",
                    high = "navy") +
  labs(title = "All Species",
       fill = "Summed \nSolution",
       x = NULL,
       y = NULL) +
  theme_minimal()

all_spec

Analyze Endangered Species (Optional)

Create Dataframes

Here we create spec and puvsp data frames similar to those we had for all species, but containing only endangered or threatened species.

end_status <- status %>% 
  mutate("endangered" = case_when(
    str_detect(status, pattern = "endangered") == TRUE ~ "yes",
    str_detect(status, pattern = "threatened") == TRUE ~ "yes",
    T ~ "no")) %>% 
  filter(endangered == "yes")

end_spec <- merge(end_status, spec, by = "id") %>% 
  select(id, amount, spf, name.x) %>% 
  rename("name" = "name.x")

puvsp_id <- puvsp %>% 
  rename("id" = "species")

end_puvsp <- merge(end_status, puvsp_id, by = "id") %>% 
  rename("species" = "id")

Run Marxan

marx_end <- marxan_problem(x = pu,
                         spec = end_spec, 
                         puvspr = end_puvsp, 
                         bound = NULL,
                         blm = 0)

marx_end_problem <- marx_end %>% 
  add_gurobi_solver(gap = 0.15) %>% 
  add_pool_portfolio(method = 2, number_solutions = 100)

marx_end_soln <- solve(marx_end_problem)

## Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (mac64)
## Optimize a model with 18 rows, 670 columns and 970 nonzeros
## Model fingerprint: 0x6e88f508
## Variable types: 0 continuous, 670 integer (670 binary)
## Coefficient statistics:
##   Matrix range     [1e+00, 1e+00]
##   Objective range  [2e+00, 5e+06]
##   Bounds range     [1e+00, 1e+00]
##   RHS range        [9e-01, 1e+02]
## Found heuristic solution: objective 4.174920e+07
## Presolve removed 14 rows and 203 columns
## Presolve time: 0.00s
## Presolved: 4 rows, 467 columns, 495 nonzeros
## Variable types: 0 continuous, 467 integer (467 binary)
## Found heuristic solution: objective 3.534249e+07
## Presolved: 4 rows, 467 columns, 495 nonzeros
## 
## 
## Root relaxation: objective 2.677119e+07, 1 iterations, 0.00 seconds
## 
##     Nodes    |    Current Node    |     Objective Bounds      |     Work
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
## *    0     0               0    2.677119e+07 2.6771e+07  0.00%     -    0s
## 
## Optimal solution found at node 0 - now completing solution pool...
## 
##     Nodes    |    Current Node    |      Pool Obj. Bounds     |     Work
##              |                    |   Worst                   |
##  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
## 
##      0     0          -    0               - 2.6771e+07      -     -    0s
##      0     0          -    0               - 2.6771e+07      -     -    0s
##      0     2          -    0               - 2.6771e+07      -     -    0s
## 
## Explored 4914 nodes (1716 simplex iterations) in 0.32 seconds
## Thread count was 1 (of 4 available processors)
## 
## Solution count 100: 2.67712e+07 2.67732e+07 2.67736e+07 ... 3.14836e+07
## 
## Optimal solution found (tolerance 1.50e-01)
## Best objective 2.677119300000e+07, best bound 2.677119300000e+07, gap 0.0000%

Sum solutions

marx_end_ssoln <- marx_end_soln %>% 
  mutate(sum = rowSums(.[6:105])) %>% 
  select(id, cost, status, locked_in, locked_out, sum)

hist(marx_end_ssoln$sum,
     main = "Selection frequencies",
     xlab = "Number of runs that the units were selected")

Visualize Results

Join polygons with output

parcels_marx_end <- inner_join(parcels, marx_end_ssoln, by = "id")

ggplot(data = parcels_marx_end) +
  geom_sf(aes(fill = sum),
          color = "white",
          size = 0.05) + 
  scale_fill_gradient(low = "slategray2",
                      high = "navy") + 
  labs(fill = "Summed \nSolution") +
  theme_minimal()

Add basemap

end_spec <- 
  ggmap(morrobay) +
    geom_sf(data = parcels_marx_end,
            aes(fill = sum),
          color = "white",
          size = 0.1,
          alpha = 0.85,
          inherit.aes = FALSE) + 
  coord_sf(crs = st_crs(4326)) +
  scale_fill_gradient(low = "slategray2",
                      high = "navy") + 
  labs(title = "Endangered & Threatened Species",
       fill = "Summed \nSolution",
       x = NULL,
       y = NULL) +
  theme_minimal()

end_spec

Combine Maps with Patchwork

This step isn’t necessary, but is really convenient if you want to present your maps together in your report. patchwork uses PEMDAS to combine ggplots together into one graphic. For example, plot_1 + plot_2 will make an image of plot_1 next to plot_2, whereas plot_1 / plot_2 will make an image of plot_1 over plot_2. For more information, see https://github.com/thomasp85/patchwork.

First, make the all_spec graph without a legend so that the combined image has only one legend.

all_no_lgnd <- all_spec +
  theme(legend.position = "none")

Next, combine the graphs with patchwork

spec_graphs <- all_no_lgnd + end_spec

spec_graphs

Finally, save the image to add to your report!

ggsave("spec-graphs.png")

Questions

The following are some questions to consider in your analysis and report.

• Are the same management units commonly included in replicate runs? You can check this by looking at the distribution of the summed solutions. What does this pattern suggest?

• Are your reserves clustered? Why or why not?

• What was the conservation goal of this preliminary analysis? Is this a reasonable goal? How may your results have changed with constraints on your conservation goal or different conservation goals altogether?

• If you did the extra analysis, how or why would sequencing of parcel protection affect future priorities?